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SECTION  I 


INTRODUCTION 

Bird  impact  hazards  to  high  speed,  low  flying,  aircraft  have  become 
one  of  the  major  flight  safety  problems  of  the  jet  age.  The  Douglas 
Aircraft  Company's  Windshield  Technology  Demonstrator  Program  is  part 
of  an  effort  by  the  Air  Force  Flight  Dynamics  Laboratory  Improved  Wind- 
shield Protection  Program  to  develop  technologies  that  will  allow  the 
design  of  aircraft  transparent  enclosures  to  increase  protection  against 
birdstrikes . 

One  of  the  basic  objectives  of  the  Windshield  Technology  Demonstrator 
Program  is  to  develop  windshield  design  technology  for  application  to 
high  performance  military  aircraft.  Pursuant  to  this  effort,  a need 
arose  for  a simplified  analytical  method  of  rapidly  assessing  the  struc- 
tural effectiveness  of  candidate  laminate  configurations. 

Typical  laminated  configurations  consist  of  alternating  structural 
plies  and  interlayers.  Structural  plies  are  composed  of  high  modulus 
materials  such  as  polycarbonate,  glass  or  acrylic.  Interlayers  are 
composed  of  materials  that  may  be  less  stiff  by  several  orders  of 
magnitude. 

A realistic  method  of  static  analysis  for  a laminated  beam  carrying 
a concentrated  load  at  the  center,  representing  a bird  impact,  and 
having  various  types  of  end  supports,  would  be  a useful  tool.  Such  a 
method  should  avoid  the  frequently  made  Bernoulli-Euler  assumption  that 
plane  sections  remain  plain  during  deformation,  because  shear  deforma- 
tions in  the  soft  interlayers  are  large  and  must  be  considered.  The 
significance  of  these  deformations  is  shown  experimentally  in  Reference 
1.  Other  rough  assumptions  that  would  unnecessarily  diminish  the  effec- 
tiveness of  the  analysis  should  also  be  avoided.  The  method  should 
accommodate  laminates  composed  of  nine  layers  or  more,  with  different 
material  properties  for  each  layer. 
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Apparently  no  method  meeting  these  requirements  exists  in  the 
literature,  although  related  work  is  described.  Reference  2 contains 
an  analysis  of  a three  layered  beam  based  on  the  assumption  that  the 
angle  of  rotation  between  a cross  section  of  the  center  layer  and  its 
undeformed  position  is  a constant  factor  times  the  slope  along  the 
length  of  the  beam.  This  assumption  is  considered  unnecessary,  and  its 
effects  on  the  reliability  of  the  results  are  hard  to  assess.  The  sol- 
ution is  based  on  the  energy  method,  which  introduces  further  approxi- 
mations. Reference  3 presents  methods  of  analysis  of  multiple  layered 
beams  based  on  the  Bernoull i-Euler  hypothesis;  consequently, the  approach 
is  not  applicable  to  beams  with  soft  interlayers.  Reference  4 offers  an 
analysis  of  three  layered  beams  involving  simplifying  assumptions  which 
are  also  not  considered  acceptable  for  the  present  application. 

Therefore, a new  approach  to  the  problem  was  developed,  and  is 
presented  in  subsequent  sections  of  this  report.  The  approach  is  based 
on  the  assumption  that  each  structural  ply  can  be  treated  as  a beam  to 
which  the  Bernoulli -Euler  hypothesis  is  applicable,  but  no  such  assump- 
tion is  applied  to  the  cross  section  as  a whole.  Structural  plies  can 
bend  and  stretch,  but  shear  deformations  of  these  plies  are  considered 
negligible.  Interlayers  are  assumed  to  carry  shear,  but  not  axial  loads 
or  bending  moments  • Deformations  through  the  thickness  of  the  beam  are 
considered  negligible,  but  stresses  normal  to  the  layers  are  assumed  to 
exist.  The  equations  of  equilibrium  and  compatibility  are  written. 

The  resulting  set  of  differential  equations  is  then  expressed  as  a single 
matric  differential  equation, which  is  solved  exactly. 

The  analytical  results  have  been  translated  into  an  efficient  Fortran 
program.  This  program  applies  to  nine  layer  laminates,  which  is  an  ade- 
quate number  in  most  cases.  The  program  can  be  easily  extended  to  cover 
more  layers,  if  necessary.  Fewer  layers  can  be  accommodated  by  intro- 
ducing negligible  stiffness  properties  (E  and  6)  for  some  of  the  layers. 
The  method  applies  to  fixed  ended  beams,  but  the  program  can  be  modified 
to  be  applicable  to  other  boundary  conditions.  Extension  of  the  method 
to  other  loading  conditions  is  possible. 
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The  most  significant  simplification  involved  in  the  analysis  is  the 
assumption  that  deformations  through  the  thickness  are  negligible. 

This  assumption  greatly  simplifies  the  analysis  without  slighting  the 
primary  feature  of  windshield  laminate  behavior,  which  is  the  relative 
freedom  of  structural  plies  to  slide  past  each  other  because  of  the  low 
stiffness  of  interlayer  materials.  The  only  negative  effect  of  the 
assumption  is  that  stresses  normal  to  the  layers  are  not  correctly  pre- 
dicted. This  is  believed  to  be  a localized  effect  confined  to  the  center 
and  ends  of  the  beam  where  loads  are  applied.  The  elimination  of  this 
assumption  is  a possible  subject  for  additional  research.  An  analysis 
which  accounts  for  transverse  deformations  might  provide  data  that  would 
be  useful  in  defining  adhesive  strength  needed  to  prevent  delamination 
in  regions  of  high  transverse  loads. 

Results  of  the  analysis  have  been  correlated  with  finite  element 
results  and  test  data  as  described  in  a subsequent  section.  The  com- 
parisons are  good  and  verify  the  validity  of  the  basic  assumptions. 

Although  the  present  method  is  intended  to  apply  to  fixed  ended 
beams,  it  can  also  be  applied  to  a beam  with  pinned  ends  by  considering 
a fixed  ended  beam  twice  as  long  as  the  beam  under  consideration.  The 
data  output  by  the  computer  program  between  the  quarter  points  of  the 
fixed  ended  beam  is  applicable  to  the  pin  ended  case. 
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SECTION  II 


ANALYSIS 


DERIVATION  OF  EQUATIONS 

This  analysis  applies  to  the  fixed  ended  laminated  beam  of  Figure  1. 
As  the  figure  shows,  L is  the  length  of  the  beam,  and  F is  the 
concentrated  load  acting  at  the  center. 

t 

I 

F 

1 

< L ► 


Figure  1.  Fixed  ended  laminated  beam. 

The  cross  section  is  shown  in  Figure  2.  The  cross  hatched  layers 
are  called  "structural  plies",  while  the  other  layers  are  "interlayers". 
The  width  of  the  beam  is  b,  and  the  thickness  of  the  ith  layer  is  t^ . 
Different  material  properties  can  be  assigned  to  each  layer.  The  flexi- 
bilities of  the  interlayers  are  assumed  to  be  large  compared  to  the 
structural  plies.  The  analysis  can  be  applied  to  a beam  having  fewer 
than  nine  layers  by  assigning  negligible  values  of  Young's  modulus,  E, 
and  the  shear  modulus,  G,  to  some  of  the  layers. 
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Figure  2.  A typical  transparency  nine  ply  cross  section. 

Because  of  symmetry,  consider  one-half  the  beam  as  shown  in  Figure  3. 
In  the  figure  i = L/2,  and  P = F/2.  The  figure  also  shows  the  xyz 
reference  frame,  and  the  definition  of  the  vertical  displacement  v. 


Figure  3.  Half -beam. 
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Assumptions 

1.  Structural  plies  carry  axial  load,  bending  moment,  shear,  and  normal 

2.  Interlayers  carry  shear  and  normal  stresses  only. 

3.  Plane  sections  of  structural  plies  remain  plane  and  normal  to  their 
elastic  axes  during  deformation. 

4.  Normal  strains  through  the  thickness  for  structural  plies  and  inter- 
layers are  negligible  and  can  be  considered  equal  to  zero. 

5.  Shear  strains  for  structural  plies  are  negligible. 

6.  All  stress-strain  relations  are  linear. 

7.  Structural  ply  bending  conforms  to  small  displacement  theory. 
Equilibrium,  Structural  Ply  Element,  i Odd 

Figure  4 shows  the  equilibrium  of  an  element  of  the  ith  structural  ply. 

, V.  and  are  the  axial  load,  shear  and  bending  moment  acting  on 

the  element  cross  section.  The  transverse  load  per  inch,  w.. , and  the 
shear  flow,  q^,  act  upon  the  lower  face.  The  corresponding  forces  on 
the  upper  face  are  w^_.|  and  q.j_-|i  since  the  i-lst  element  is  next 
above.  Note  that  i is  odd  for  structural  plies. 
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load. 


I 


1 

1 


q.dx 


w^dx 

■< dx  ► 

Figure  4.  Structural  ply  equilibrium. 

Summing  forces  in  two  directions  and  moments  about  0 gives: 


dx 

+ 

'’i-1 

- ‘’i  = 

0 

(i  odd) 

dV. 

dx 

+ 

^•-1 

- ”i  = 

0 

(i  odd) 

dM. 

dx 

+ 

- 

Vi^ 

- 

^i 

Sr  = 0 

(i  odd) 

0) 

(2) 

(3) 


where  1/2  dV^  has  been  deleted  from  Equation  3 to  reduce  complexity 
of  problem.  This  is  allowable  due  to  insignificance  of  the  term. 
Equilibrium,  Interlayer  Element,  i Even 


Figure  5 shows  the  equilibrium  of  an  element  of  the  ith  interlayer. 
The  notation  is  consistent  with  the  notation  established  for  structural 
plies,  although  no  axial  loads  or  bending  moments  are  acting,  because 
of  the  assumption  that  these  forces  are  negligible  for  interlayers. 

Note  that  i is  even  for  these  plies. 
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Figure  5.  Interlayer  equilibrium. 


Sunwing  forces  in  two  directions  and  moments  about  0 gives; 


o 

M 

O' 

1 

1 

cr 

(i  even) 

(4) 

1 '“''i 

(i  even) 

(5) 

( r ■ '’i  ° 

(i  even) 

(6) 

From  Equation  4,  = 'll  » ^^liJation  6,  V.  = q.t. 

• 

Therefore,  from  Equation  5, 

dq. 

t.  w.  , - w.  = 0 

1 dx  1-1  1 

(i  even)  . 

(6a) 

Also  note  that 

qo  = qg  = 0 

(6b) 

o 

II 

II 

o 

(6c) 

since  transverse  forces  and  shear  flows  are  assumed  to  be  zero  on  the 
upper  and  lower  surfaces  of  the  beam. 
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Compatibility,  Structural  Plies,  i Odd 

Figure  6 shows  the  relationships  between  the  displacements  of  the 
centerline  of  a structural  ply  and  the  displacements  of  the  upper  and 
lower  surfaces  of  the  same  layer.  In  the  figure,  u^  , and 
are  the  longitudinal  and  transverse  displacements  and  the  slope  of  the 
centerline  of  the  ith  ply.  The  longitudinal  and  transverse  displacements 
of  the  lower  surface  of  the  ith  ply  are  denoted  and  . Conse- 
quently the  corresponding  displacements  of  the  upper  surface  are  n^_i 
and  ^ , since  the  upper  surface  of  the  ith  ply  is  the  lower  surface 
of  the  i-lst  ply. 

Displaced  ply  axis 


Figure  6.  Compatibility. 

t^. 

Then  from  Figure  6,  ^i-1  " ^i  " ®i 

t. 

and  + e^.  ^ for  small  di splacements , (8) 

These  equations  are  consistent  with  the  assumption  that  cross  sections  of 
the  structural  ply  remain  normal  to  tne  elastic  axis  during  defonnation. 
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Also  from  Figure  6, 


(9) 


i-1 

= "i  ■ 

^i 

= ''i  ’ 

dVi 

®i 

(10) 

(11) 


Equations  9 and  10  hold  because  strains  through  the  thickness  of  structural 
plies  are  considered  negligible.  Equation  11  is  consistent  with  the  assump- 
tion that  structural  ply  bending  conforms  to  the  theory  of  small  displacements. 


Compatibility,  Interlayers,  i Even. 


Strains  through  the  thickness  of  the  interlayer  are  assumed  negligible. 


i-1 

= • 

02) 

In  consequence  of  Equations  9, 

10 

and  12: 

to  = vi  = Cl  = = ''3  = 

^3 

* C4  = ••• 

= ^7 

fl 

CD 

II 

< 

VO 

It 

^9  • 

V.  - V 

i 

= 1,  3,  ... 

9 

(i  odd) 

03) 

ti  = V 

i 

* 0,  1.  2 . 

..  9 

(all  i) 

04) 

e.  = = e 

^ dx 

i 

= 1,3  .. 

. 9 

( i odd ) 

05) 

All  vertical  displacements  and 

slopes  are  the 

same. 

Strains,  Structural  Plies,  i 

Odd 

• 

Longitudinal 

displacement: 


strain,  e,  is  equal 


to  the  derivative  of  longitudinal 
axial  strain  . 
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In  small  displacement  theory,  curvature  B is  equal  to  the  rate 
of  change  of  slope,  therefore. 


The  notation,  dlT  ~ ^ ’ 05b) 

to  denote  the  rate  of  change  of  beam  curvature,  is  useful  subsequently. 


Strains,  Interlayers, 


Figure  7 shows  the  relationships  between  the  shear  strain  of  the  ith 
interlayer,  , the  slope  of  the  beam,  6 , and  the  longitudinal 

displacement  of  the  ith  interlayer  upper  and  lower  surfaces.  In  the 
figure,  denotes  the  rotation  of  a cross  section  of  the  ith  inter- 
layer relative  to  its  undisplaced  position. 


Figure  7.  Interlayer  shear  strain. 


Force-Deformation  Relations 


The  longitudinal  force  on  the  ith  interlayer  is  given  in  terms  of 
the  longitudinal  strain  by 

H,  = E,  A,  e . 


H,  ' E,  A 


du  • 


i i dx 


(i  odd) 


(16) 


where  and  are  the  values  of  Young's  modulus  and  the  cross- 
sectional  area  of  the  ith  structural  ply.  The  bending  moment  in  the 
ith  structural  ply  is  given  in  terms  of  the  curvature  by 


Ei  I.  6,. 


d^v 


(i  odd) 


(17) 


where  is  the  moment  of  inertia  of  the  cross  section  of  the  ith 

structural  ply.  The  shear  flow  in  the  ith  interlayer  is  given  in  terms 
of  the  shear  strain  by 


G.  b Yi  . 


where  is  the  shear  modulus  of  the  ith  interlayer. 

Matric  Formulation 

Equation  1 can  be  written  for  each  structural  ply  thus; 

dH, 

dT  * '<0  ■ ‘'1  ' “ 

dH. 
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Equation  18a  has  been  derived  in  detail  to  show  the  method  of 
converting  the  scalar  differential  equations  into  matric  form.  Subse- 
quent matric  equations  are  derived  in  a similar  manner,  although  the 
details  are  omitted. 


The  matrices  in  Equation  18a  are  denoted  according  to  the  symbols 
written  under  the  equation.  Thus,  H is  a column  matrix  of  structural 
ply  longitudinal  forces,  dH/dx  is  its  derivative,  q is  a column 
matrix  of  interlayer  shear  flows,  and  is  a Boolean  differencing 
matrix. 

Therefore , 

^ % 1 ■ 0 . t'9) 


From  Equations  3,  4 and  6b, 


The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  M denotes  a column  matrix  of  structural  ply 
bending  moments,  dM/dx  is  the  derivative  of  M,  V is  a column  matrix 
of  structural  ply  shear  forces,  and  is  a rectangular  matrix  of 

structural  ply  thicknesses. 


V 


0 . 


(22) 
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From  Equations  2,  6a  and  6c: 


The  matric  fonti  of  this  equation  is: 


1 1 1 M 


dV-j/dx 
dV3/dx 
dVg/dx 
dVy/dx 
.dVg/dx 

dx 


dq2/dx 
dq^/dx 
j dqg/dx 
(^dqg/dx 


'I 


The  matrices  are  denoted  according  to  the  symbols  written  under  the  equa- 
tion. The  matrices  dV/dx  and  dq/dx  are  the  derivatives  of  the  column 
matrices  V and  q previously  defined,  C is  a Boolean  column  matrix, 
and  tj  is  a column  matrix  of  interlayer  thicknesses.  The  superscript  T 
indicates  a transposed  matrix. 

-tT  dV  . 


dx 


dx 


(23a) 


From  Equation  16: 


/ ^1  \ ” 

-EiA^ 

''du^/dx\ 

1 ^3  1 

E A 

dug/dx  1 

E5A5 

< 

dug/dx  \ 

/ ^7  ( 

E7A7 

du^/dx  1 

\ 9 / 

EgAg^ 

^dug/dx  J 

du 


The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  u denotes  a column  matrix  of  longitudinal  dis- 
placements of  structural  ply  centerlines,  du/dx  is  the  derivative  of 
u,  and  is  a diagonal  matrix  of  structural  ply  axial  stiffnesses. 

• H = k — (24) 
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From  Equation  17  : 


E I \ 

/A  ^ 

hh 
hh 

E9I9 
M ki 

The  matrices  are  denoted  according  to  the  symbols  written  under  the  equa- 
tion. The  symbol  kj  denotes  a column  matrix  of  structural  ply  bending 
stiffnesses . 

d^v 


M 


^1  —7 
dx^ 


(25) 


From  Equation  18: 
‘^2 


82b 

G^b 

Ggb 


dv 
1 3)T 


r ^ 


G4b 


Geb 


Gsb 

^8  J 


’^1-^2 

Inj-n^ 

n^-ng 


An 


The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  ^ denotes  a column  matrix  of  longitudinal 
displacement  differences  at  layer  interfaces,  k^^^  is  a column  matrix 
of  interlayer  shear  stiffnesses,  and  k^  is  a diagonal  matrix  of 
interlayer  shear  stiffnesses  divided  by  interlayer  thicknesses. 


q 


k An 
n 


But 

where 


An  = A n 
n 


1 -1 


1 -1 


1 -1 


1 -1 
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and 


n = I I ri2  I TI3  I I ng  I Tig  I ny  I ng  I (column)  . 

An  is  a rectangular  differencing  matrix,  and  n is  a column 
matrix  of  longitudinal  displacements  of  layer  interfaces. 


k ^ + k An 
t dx  n 


The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  n denotes  a column  matrix  of  longitudinal  dis- 
placements, t^  IS  a column  matrix  of  structural  ply  thicknesses,  and 
and  are  Boolean  coefficient  matrices. 
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From  Equation  24 


0 


(28) 


where 

Eliminate 

But  k^= 
where 


and  tj  is 
where 


El iminate 

where 

and 


^ - a,  , H = 
dx  1 ,2 

®1,2  ^ ’'a  • 

n from  Equations  26  and  27: 

‘n^t  ’ " \ • ^^u  = ‘^5  . 

^(t^  + t3) 

^(ts  + t^) 

7(t7  + tg) 

1 column  matrix  of  average  structural  ply  thicknesses. 


(29) 


Z^ti+t2+jt3' 

ita+t^+^sj 


2t7+ 


1, 


(30) 


q and  ^ from  Equations  15,  19  and  30 

S'  ^2,1  ^ ■ ^2,3  ® ^ 


a,  , = A-  k A„ 

2,1  q n q 


32,3 


-A  k t 

q n 


(31) 

(32) 

(33) 
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Eliminate  V from  Equations  22  and  23a. 


(34) 


since  = t^.  Eliminate  M from  Equations  25  and  34, 


^ - '‘l  & ■ 0 


where  K 


J = (EI)^  + {ED3  + (EDg  + (ED7  + (El)g 


(35) 

(36) 


Therefore,  from  Equations  15,  15a,  15b,  and  35 

S - "i  ^ = “ ■ 


(37) 


Eliminate  q from  Equations  30  and  37.  Then, 


^ . K-''  ^ - a ^ = 0 

dx  ®2,3  dx  ^5,4  ^ 


a.  . = k:^  t"'^  k t 
5,4  I n 


where 

Therefore,  from  Equations  15,  15a,  28  and  38, 


where 


dx  ■ ®5,2^  ‘ ®5,4  6=0* 
®5,2  " '^I^®2,3®1,2  ' 


(38) 
(38a) 

(39) 

(40) 


Equations  15,  15a,  I5b,  28,  31  and  39  are  a set  of  14  linear,  ordinary, 
first  order,  homogeneous  differential  equations  with  constant  coeffi- 
cients, in  the  14  dependent  variables  u^,  H^,  v,  e , 6 and  ^ . 

The  number  of  variables  is  14  since  the  matrices  u and  H each  con- 
tain five  variables.  The  theory  of  differential  equations  shows 
that,  for  a set  of  equations  of  this  kind,  n linearly  independent 
solutions  exist,  where  n is  the  number  of  equations,  equal  to  14 
in  this  case.  See  Reference  5,  article  12-8. 


20 


The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  Y denotes  a column  matrix  of  beam  responses, 
dY/dx  is  its  derivative  with  respect  to  x , and  A is  a square 
matrix  of  constant  coefficients. 

^ - ay  = 0 (42) 

This  linear  matric  differential  equation  can  be  solved  by  taking 

Y = Ge^’^  (43) 

where  G is  an  unknown  column  matrix,  X is  an  unknown  scalar  and 
e is  the  base  of  natural  logarithms.  Eliminating  Y from  Equations 
42  and  43  gives 

(A-XI)G  = 0 . (44) 

Evidently,  x and  G are  an  eigenvalue  and  an  eigenvector  of  the 
characteristic  Equation  44.  Therefore  the  solution  of  Equation  42 
can  be  written 

14  \.x 

Y = E C.  G.  e " (45) 

k=l  ^ 

where  Gj^,  and  are  the  kth  eigenvector  and  eigenvalue  of  Equation 
44,  and  the  Cj^'s  are  arbitrary  constants  needed  to  satisfy  the  boundary 
conditions. 


Nature  of  the  Roots  of  the  Characteristic  Equation 


The  eigenvalues  (roots)  of  Equation  44  are  expected  to  be  all  real, 
from  the  nature  of  the  physical  problem.  The  roots  can  be  easily  shown 
to  occur  in  pairs  equal  in  magnitude  and  opposite  in  sign,  such  that 
if  and  are  a pair  of  roots,  then 

Appendix  A.)  It  is  shown  subsequently  that  at  least  six  repeated  roots, 
equal  to  zero,  exist.  No  more  than  six  such  roots  are  expected, 
because  of  the  nature  of  the  problem. 


Algebraic  Modes 

The  solutions  of  Equation  42  of  the  form  given  by  Equation  43  are 
here  called  exponential  solutions  or  exponential  modes.  The  following 
paragraphs  show  that  other  solutions  exist  which  are  algebraic  functions 
of  X . Consequently,  these  solutions  are  called  algebraic  solutions 
or  algebraic  modes.  Six  linearly  independent  algebraic  modes  are  shown  to 
exist.  These  solutions  are  also  linearly  independent  of  the  exponen- 
tial solutions,  since  they  are  algebraic.  Therefore,  the  character- 
istic equation  (Equation  44)  can  yield  at  most  only  eight  linearly 
independent  exponential  modes.  Consequently,  six  of  the  roots  of 
Equation  44  must  be  zero. 

In  view  of  the  existence  of  the  algebraic  modes,  the  complete 
solution  of  Equation  42  can  be  written  in  the  form 

8 XkX  6 

Y * S C G e <<  + X Z,  Y„  (46) 

k=l  ^ k=l  K 

where  the  C/\  's  are  arbitrary  constants,  and  the  Ya  's  (algebraic  modes) 

are  column  matrices  whose  elements  are  algebraic  functions  of  x.  The 

only  requirements  for  the  Y/\  's  is  that  they  must  satisfy  Equation  42, 

and  the  G^'s  and  Ya  's  must  be  linearly  independent.  The  following 
K K 

Yft  's  satisfy  these  requirements: 

In 
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Rigid  Body  Translation  Parallel  to  x - 


where  u 


5 rows 


(47) 


Uniform  Extension  Parallel  to  x 

% V(o)i  M “ r I 


Rigid  Body  Translation  Parallel  to  z 


5 5 

rows  rows 


0 I 0 I o|  (Column) 


(49) 


Rigid  Body  Rotation  about  y 


= {u(T)  I 0 I X I 1 I 0 I o}  (Column)  (50) 
5 rows 


where  u^^  = Qt  + ®^(0) 


-l-l-l-l 

-1-1-1 

-1-1 

-1 

0 0 0 0 


a = a scalar 


(51) 


Bending  about  y - 

V I"  r I “ ^ 
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•*» 


Shear  Parallel  to  z - 


V " { “(2)  V(l)’^  6 2 


I (Column) 


where  ^ column  matrix  with  5 rows, 

o and  U/»x  must  satisfy  the  equation 


Su^2)  - Ta  = W 
where  S = . T = k^u^^j  . W = k^Qt 


The  matrix  S is  singular.  Equation  54  is  a matrix  equation 
equivalent  to  five  scalar  simultaneous  equations  in  the  five  unknown 
elements  of  u^2)  unknown  o.  However  one  element  of  u^2) 

be  chosen  arbitrarily  since  S has  one  degree  of  singularity.  In 
this  respect  u^2)  character  of  an  eigenvector.  Partition  S, 

T,  W and  U/o\: 


^ ' ^1 ,1  ^1 ,2  ^ ^ ^ 
^2.1  ^2.2 


^ ) 1 row 
T2  ] 4 rows 


W = ' ) 1 row 

, W2  J ] 4 rows 


1 C0I . 4 col. 


“(2)  " I ° |M  <"Ow 
I U(2)  ) } 4 <"Ows 


where  the  element  in  the  first  row  of  t>een  chosen  equal  to  zero. 

From  Equations  54,  56  and  57, 

^1,2  ^(2)  ■ ’’’f  " *^1 
■ ^2,2  ‘^(2)  "V  " '^2 
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a 


(59) 


" Tt  [^1.2  ^2)  ■ ‘^l] 
and  U(2)  = (52^2  " ^2  ^1,2^  ^ ^^2  ' ^2 

Nature  of  the  Algebraic  Solutions 

The  first  five  algebraic  solutions  (Y/\  , Y/\  , ...  Y/\  ) correspond 

• 2 ^ 

to  solutions  provided  by  the  engineering  theory  of  tension  and  bending 
of  beams.  These  solutions  have  been  appropriately  named  to  reflect 
this  fact.  Thus  Y;\^  can  be  recognized  as  corresponding  to  a rigid  body 
translation  of  the  beam  parallel  to  x,  since  the  longitudinal  motions 
of  the  centerlines  of  the  structural  plies  are  all  equal  to  1,  and  all 
other  responses  are  zero.  The  mode  Y^^  can  be  recognized  as  a rigid 
body  rotation  about  the  centerline  of  the  bottom  structural  ply,  repre- 
sented by  Qt  , plus  a rigid  body  translation  parallel  to  x , repre- 
sented by  The  mode  Y/^^  corresponds  to  the  responses  calculated 

by  engineering  theory  when  the  beam  is  subjected  to  a constant  shear 
load.  However  the  responses  given  by  elementary  theory  involve  an 
approximation,  whereas  the  responses  given  by  exactly  satisfy 

the  differential  equations. 


Response 

Appendix  A shows  that  the  eigenvalues  occur  in  pairs,  such 

that  for  a pair  x and  x.  , 
a 0 


X 


6 


-X 


(61) 


Appendix  A also  shows  that  the  corresponding  eigenvectors  have  the 
property  that  if 
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(62) 


where  G , Gu  etc.  are  partitions  of  G„  corresponding  to  the 

Uq  “a  ^ 

partitions  u,  H etc.  of  the  matrix  Y (Equation  41).  Because  half  of 

Xx 

the  eigenvalues  X are  positive,  the  corresponding  functions  e 
can  become  very  large  as  x increases.  To  avoid  this  difficulty, 
proceed  as  follows:  See  Equation  46.  Let 


8 X.x 

Y,  = Z C,  e ' 


k=l 


(exponential  response) 


= Y + Y 
^Et  '^Eg 

where  4 

Yr.  = E G„  e , 
k=l  'k  \ 


X^kX 


(63) 

(64) 


and  where  x^j^  and  Xg^^  are  a pair  of  eigenvalues  having  the  properties 


^k  ^k 

vectors.  Now 

9k 

0 . 

G.  and  G 
^ 9k 

are  corresponding  eigen- 

''Eg 

II 

9k  9k 

Xg  (X-t+t) 

e 

(65) 

II 

Z G 
9k  9k 

*^g|^(t-x)^Xg|^S, 

(66) 

II 

9k  9k 

^ tt ( ^”X  ) 

e 

(67) 

where  Cg|^ 

" ^9k 

Xgkt 

e 

(68) 
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From  Equation  46: 


4 Xji  X 4 6 

V = £ \ * E,  e " * E c,  y,,  . (69) 


k=l  k 


^k  k=1  "k  "k 


Thus  the  troublesome  eigenvalues  , which  are  greater  than  zero, 

have  been  lumped  into  the  arbitrary  Constants  Cn  , and  thus  eliminated. 
^XlX  Xj,.  (fc-x)  k 

The  functions  e and  e are  less  than  or  equal  to  1 for 

0 < X < t , since  X..  < 0.  Let 


^k 

^tbX  Xj,i  (t-x) 

F (x)  = e and  F (x)  = e ^ 

*k  9k 


(70) 


4 4 6 

••  » = E G.  F (x)  C + E G F (x)  c * E C,  y,  (71) 

k=l  ‘k  'k  ‘k  k-1  Gk  9k  S,  fl,  \ 


y(x)  . H,  F^„(x)  C,  + Hg  Fg„(x)  Cg  k y^(x) 


where  the  symbols  in 
G. 


= Fg. 

I 


(72) 

this  equa^on  are  defijied  by  Equations  73^74,75  and  7 

(73) 


3 4 


= C 

ij  9 [_ 


G G 
9l  92 


% \ 


S.D' 


(74) 


/ c„  \ C = I 

f C ^ 

c,  = 

c.  \ 

( M 

9l 

A 

’ 1 

1 0 1 

1 c 

c.  / 

1 ^2  ( ) 

92 

> 

^2  ( 

c 1 

• / 

1 ^ 1 

93  i 

• \ 
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V ^4  / 

^ 94  ^ 
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The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  partitions  of  Y;\(x)  are  defined,  from  Equation  76,  by  the 
following  equations; 


“(0) 

“(0)"  “ “(1)  "(D* 

AH<”>  ' 

V(o)  ° ° V(i)  Vo)”  1 

Av<«>  • 

[“ 

0 ^ 1 1 ] 

A«'*>  ■ 1 

C\J 

CM 

X 

X 

o 

o 

A8<«>  = 

0 

o 

o 

o 

X 

A,(«>  ■ 

0 

0 0 0 0 1 J 

(78) 

(79) 

(80) 
(81) 
(82) 
(83) 


Interlayer  Shear  Flows,  Structural  Ply  Shears  and  Moments 


From  Equations  15  and  30 

q = k^te  - u . (84) 

From  Equations  15b,  22  and  25, 

V = V^q  - (85) 

where  . 

From  Equations  15a  and  25 

M = kj6  . (86) 

Let  S(x)  = total  shear  on  the  cross  section  at  station  x. 


The  matrices  are  denoted  according  to  the  symbols  written  under  the  equation. 
Therefore  S = u|o) 

Eliminate  V from  Equations  85  and  88, 

S = t^  q - Kj(^  . (89) 

Eliminate  q from  Equations  84  and  89, 

S * a2  3 u + Kj  a^  ^ e - Kj4i  . (90) 
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Figure  8.  Boundary  conditions. 

For  all  X - 

The  total  shear  on  the  cross  section  is  constant 

S(x)  = P . 

S(x)  is  most  easily  evaluated  at  X = 0,  therefore  take 
S(0)  = P (1  condition)  . 

From  Equation  90, 

-Kj^(O)  = P and  «(0)  = -k‘V  , 

since  u(0)  = 6(0)  = 0. 

Condition  1.  Fixed  Ends  at  x = . For  x = *•  — 

j u(t)  =0  (5  conditions) 

I 6(t)  =0  (1  condition) 

1 


(92) 


(93) 


From  Equation  77 


* "...f"n(0)C„  ♦ V,„(0)C 


ut  gD' 


- «vtf90"»S  * ''a»<">'a  = 
* "e/gD«»S  * W°>':a 

"u/tO<‘>':£  * "u/gD">S  * ''au*'>'^A  = 
* "etfgD<'>S  * ''a6<'>':a  = 


-k:'p 


"u.fgo"” 

^U 

v/aO<0) 

-"vA^goW) 

%<V'» 

\.'gD<“> 

-e- 

<C 

>- 

hu 

%/gO<A> 

he 

The  matrices  are  denoted  according  to  the  symbols  written  under  the 
equation.  The  symbol  C denotes  a column  matrix  of  arbitrary  constants, 
ij;  is  a column  matrix  of  boundary  conditions,  and  B is  a coefficient 
matrix. 

BC  = i|;  . (96) 

From  Equations  7C  and  74 

^D«»  = 

^d(')  ■ fggO'  <«8> 


r 

I 


where  I is  a unit  (identity)  matrix.  Therefore;  B = 


“(0) 

0 

0 

‘'(I) 

0 

'"(Z) 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

“(0) 

‘"(0)^ 

0 

‘'(I) 

U(i)t 

0 

0 

0 

1 

1 

2 

...  (99) 


Effective  Beam  Stiffness 

For  a fixed-ended  monolithic  beam  carrying  a concentrated  load 
at  the  center, 


v(a)  - ^2EI  • 


(TOO) 


Therefore  define  the  effective  stiffness  of  the  laminated  beam  as 


P!.3 

^^eff  ' I2v(t) 


(101) 
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EQUATION  SUMMARY 


This  section  surmiarizes  the  equations  for  the  laminated  beam. 
Figure  8a  shows  the  basic  beam  dimensions. 


F/2  F/2 


Side  View 


Figure  8a.  Beam  dimensions. 

Odd  numbered  layers  are  structural  plies. 

Even  numbered  layers  are  interlayers. 


Input 


F,  I , b , X = X.  , j = 0,l,2...n  where  x = 0 , x 
J ^ 0 fi 


For  the  structural  plies(i  odd); 


For  the  interlayers  (i  even): 


ti  . E. 


ti  . G. 


Equations 


For  i odd:  A^.  = bt^. 


T = ^^^i 

1 IT 

G»b 


^ih  ^3^3  ^5^5  hh  * ^9^9 


^5‘'^64^7 


2^7^V2^9  ' 


1 -1 
1 


^1.2  ' '^a' 

^2,1  ' svj 

a,  - = -A„k  t 

2.3  q n 

®5,2  " ’^1^4,3®!  ,2 

a,.  . = K"^t^k  t 

5.4  i ri 


A = 
14x14 


0 

^1.2 
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0 
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0 

0» 

ro 

A# 

CO 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

®5.2 

0 

0 

®5.4 

0 

(106) 


(107) 


Find  the  eigenvalues  and  eigenvectors  of 


(A  - XI)G  = 0 


(108) 


Six  of  the  eigenvalues  should  be  zero.  The  remaining  eight  eigenvalues 

should  be  real  and  should  occur  in  pairs  such  that  for  any  pair, 

X.  and  X„  , 
t g ’ 

^»<0  , X >0  and  ~ ^ • (109) 

If  the  eigenvalues  do  not  occur  in  this  manner,  print  an  error  message. 
Otherwise,  discard  the  eigenvalues  which  are  equal  to  or  greater  than 
zero,  and  the  corresponding  eigenvectors.  Denote  the  remaining  eigen- 
values and  eigenvectors  X and  G , i = 1,2, 3, 4. 

*i  *"i 
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FT  ~ ■ 

Form 


^D<‘>  ■ 


S ^ 

Partition 
S = 

Discard 


(110) 


diagonalized.  Partition 


ut 
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Form  u 


q n q 
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^ “(0)  ’ 


Q = 


-1  -1  -1  -1 

-1  -1  -1 

-1  -1 
-1 

0 0 0 0 


(111) 


w = Q t . 


(112) 


S,  T and  W : 


n.l 


'2.1 


'1.2 


'2,2 


1 row , T = ( T,  ^ 1 row , W = Wl  row 


4 rows 


4 rows 


1 col . 4 eels. 


.1  ^2.1  • 


"(2) 


" ^^2,2  ■ ^1,2^  ^ ^*^2  ■ ''^2  *^1  ^ 


(114) 
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Then 


'ft/Xj)' 

fglXj)  = 

(123) 
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F,(Xj)  and  Fg(x.) 


diagonalized. 
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Partition  C : 


C 


4 rows 
4 rows 
6 rows 


'n(x.)' 
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i “("j'  ■ 

H(Xj)  , 

v(x_j) 

, e(Xj)  , 6(Xj)  , 4i(Xj) 

q(Sj)  . 

V(xj)  , 

H(xj) 

for  j =0,  1 , 2 ... 

(E')eff 

ILLUSTRATIVE  EXAMPLES 

Four  sample  cases  have  been  selected  to  illustrate  the  usefulness  of 
this  program.  Cross  sections  and  material  properties  for  the  four  cases 
are  shown  in  Figures  9,  10,  11  and  12.  The  end  conditions  have  been 
defined  as  fixed.  The  load  is  one  pound  and  has  been  applied  at  the 

beam  center  line.  The  outer  glass  shield  is  fixed  at  the  ends,  although 

for  a typical  windshield  the  glass  floats  on  an  interlayer  and  is  free- 
floating  at  the  ends.  The  program  is  set  up  to  handle  a total  of  nine 

plies.  For  beams  with  less  than  nine  plies,  such  as  Cases  B,  C and  0, 

small  values  were  assigned  to  the  thickness,  modulus  of  elasticity,  and 
shear  modulus  to  negate  the  influence  of  those  fictitious  plies. 

Input  data  and  output  data  for  five  examples  are  presented  in 
Appendix  D.  Case  A was  analyzed  for  a beam  ten  inches  long  and  for  a 
beam  30  inches  long.  Case  D and  Case  C were  analyzed  for  a beam  ten 
inches  long.  Case  B was  analyzed  for  a beam  20  inches  long.  Table  1 is 
a summary  and  explanation  of  input  data  for  Case  A.  Table  2 is  a summary 
and  explanation  of  output  data  for  Case  A.  Table  3 lists  a summary  of 
deflections  for  the  five  examples  as  well  as  additional  deflections  for 
the  four  cases  at  various  lengths  measured  at  the  beam  center  line. 

The  deflections  do  not  vary  with  the  cube  of  the  length  and  this  devia- 
tion is  shown  in  Figure  13,  which  compares  a laminated  beam  to  a multi- 
ply beam  without  interlayer  for  Cases  A,  B,  C and  D. 

At  the  end  of  output  data  for  each  example  in  Appendix  Dt  there  is  a 
printout  for  EI(EFF)  which  is  defined  as  "effective  stiffness".  This 
value  is  calculated  from  the  deflection  at  the  beam  center  line  using 
the  formula  for  a beam  with  fixed  ends  and  a center  point  load  which 
is: 


r 


END  COND  FIXED 
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CASE  A 


PLY  NO. 

1— 

E 

G 

TYPE 

1 

0.100 

10  X 10^ 

4 X 10^ 

GLASS 

2 

0.250 

183 

73 

CIP  SILICONE 

3 

0.250 

340,000 

109,000 

POLYCARBONATE 

4 

0.100 

183 

73 

CIP  SILICONE 

5 

0.250 

340,000 

109,000 

POLYCARBONATE 

6 

0.100 

183 

73 

CIP  SILICONE 

7 

0.250 

340,000 

109,000 

POLYCARBONATE 

8 

0.100 

183 

73 

CIP  SILICONE 

9 

0.250 

340,000 

109,000 

POLYCARBONATE 

SECTION  A-A 


X(0)  = 0 

DEFLECTED  SHAPE  OF  BEAM  - NO  SCALE 
See  Tables  1 and  2 for  symbols 

Figure  9.  Cross  sectional  properties  for  Case  A 


X(0)  = 0 

DEFLECTED  SHAPE  OF  BEAM  - NO  SCALE 
See  Tables  1 and  2 for  symbols 

Figure  11.  Cross  sectional  properties  for  Case  C. 
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COMPUTER  PROGRAM 
SYMBOLS 


DESCRIPTION 


INPUT  DATA 


NX 

Number  of  Increments 

10.0 

F 

Load  (Pounds) 

1.0 

L 

Half  Length  (Inches) 

5.0 

BB 

Width  (Inches) 

1.0 

T 

Ply  Thickness  (Inches) 
(5  Structural  Plies  and 
4 Interlayer  Plies) 

.10  .25,  .25,  .10,  .25, 
.10,  .25,  .10,  .25 

E 

Modulus  of  Elasticity 
(PSI) 

10^,  183,  340,000,  183, 
340,000,  183,  340,000, 
183,  340,000 

G 

Shear  Modulus  (PSI) 

4 X 10^,  73,  109,000,  73, 
109,000,  73,  109,000,  73, 
109,000 

X 



Distance  from  Center  of 
Beam  (Inches) 

0,  .50,  1.0,  1.5,  2.0, 
2.5,  3.0,  3.5,  4.0,  4.5, 
5.0 

TABLE  2.  SUMMARY  OF  OUTPUT  DATA  FOR  CASE  A 


COMPUTER  PROGRAM 
SYMBOLS 

DESCRIPTION 

OUTPUT  DATA 

X(0) 

Increment  Number 

0 

U 

Horizontal  Displacement  of 
Center  of  Structural  Plies 
( Inches) 

0 

H 

Axial  Load  in  Structural 
Plies  (Pounds) 

-.0959,  -0979,  -00192, 
-.00389,  -.1918 

VEE 

Vertical  Deflection 
(Inches)  (Zero  Deflection  is 
Point  of  Maximum  Deflection 
at  Center) 

0 

QXJ 

Shear  Flow  in  Interlayer 
(Pounds  Per  Inch) 

0 

VXJ 

Shear  Force  in  Ply 
(Pounds) 

.16,  .085,  .085,  .085, 
.085 

MXJ 

Moment  on  Structural  Ply 
( Inch-Pounds) 

.322,  .171,  .171,  .171, 
.171 

THETA 

Slope  (Radians) 

-5.551115E-17 

BETA 

Rate  of  Change  of  Slope 

.000386 

PHI 

Rate  of  Change  of  BETA 

-.000192 
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TABLE  3.  DEFLECTION  - (INCHES) 


OUTPUT 

FROM  PROGRAM 

LENGTH  OF  BEAM  - INCHES 

10 

20 

30 

50 

CASE  A 

4 Ply  Polycarbonate 

.00153 

.00733 

.0152 

.0339 

CASE  B 

2-Ply  Glass 

.0000239 

.000175 

.000520 

.00178 

CASE  C 

1-Ply  Polycarbonate 

.000252 

.00178 

.00511 

.0172 

CASE  D 

2-Ply  Polycarbonate 

.000364 

.00142 

.00292 

.00731 
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Figure  13.  Deflection  versus  beam  length 


where  F is  load,  L is  length  and  V is  deflection.  Effective  stiff- 
ness varies  with  length.  Table  4 is  a comparison  of  effective  stiffness, 
at  various  lengths  for  the  four  cases,  to  the  stiffness  if  all  the  plies 
acted  as  individual  beams.  The  reason  why  effective  stiffness  increases 
with  length  is  that  the  interlayers  are  more  effective  in  transferring 
shear  between  structural  plies,  as  length  increases.  Thus,  the  behavior 
of  the  longer  beams  approaches  monolithic. 

Figure  14  is  a plot  of  strain  (inches  per  inch)  versus  thickness 
(inches)  for  Case  A (4  ply  polycarbonate)  at  the  center  of  the  beam 
with  a one  pound  load  applied  at  the  center  (beam  length  is  ten  inches). 

The  values  were  calculated  from  the  output  data  using  the  moment  in  the 
structural  plies  (MXJ)  and  the  axial  load  in  the  structural  plies  (H). 

Figure  15  shows  a finite  element  model  for  a beam  ten  inches  long 
with  a cross  section  identical  to  Case  A.  The  beam  was  analyzed  using 
the  computer  aided  Structural  Design  computer  program*to  provide  a com- 
parison. The  deflections  agreed  within  one  percent.  This  comparison  is 
presented  as  verification  of  the  validity  of  the  basic  assumptions  for 
the  enclosed  derivations. 

An  actual  test  (Reference  1)  was  conducted  to  determine  the  deflection 
of  the  midpoint  of  a beam  with  a cross  section  similar  to  Case  C,  Figure  11. 
The  beam  length  was  34,7  inches.  The  face-ply  (ply  No.  1)  was  "free 
floating".  The  deflection  of  the  test  beam  was  .00882  inch.  The  data 
output  for  the  program  is  given  in  Appendix  0 and  shows  a deflection  of 
0.00913  inch  which  agrees  within  3.5  percent.  This  comparison  is  pre- 
sented as  additional  verification  of  the  validity  of  the  basic  assumptions 
for  the  enclosed  derivations. 


TABLE  4.  EFFECTIVE  STIFFNESS  (POUNDS-INCHES  SQUARE) 


OUTPUT  FROM  PROGRAM 

SINGLE  ACTING 
BEAMS 

(NO  INTERLAYER) 

-n 

LENGTH  OF  BEAM  - (INCHES) 

10 

20 

30 

. 

50 

CASE  A 

3,397 

5,683 

9,238 

19,212 

2,600 

4-Ply  Polycarbonate 

CASE  B 

217,500 

237,900 

270,600 

366,000 

209,100 

2-Ply  Glass 

CASE  C 

20,700 

23,450 

27,500 

37,900 

19,720 

1-Ply  Polycarbonate 

CASE  D 

14,326 

29,370 

48,175 

r ' 

89,080 

7,900 

2-Ply  Polycarbonate 

1 


« ■ 


LENGTH  = 10.0  INCH  • 

LENGTH  = 50.0  INCH • 


-31.3  X 10"^  -19.4  X 10'^ 


Figure  14.  Strain  (In/In)  versus  thickness  for 
Case  A at  center  of  beam. 
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.050 


. i^4  4 PLCS 
.053  4 PLCS 


DEFLECTION  AT  MIDDLE  OF 
BEAM  = .00153  INCHES 


Figure  15.  Finite  element  model  of  a fixed  ended  beam  - Case  A. 


SECTION  III 

COMPUTER  PROGRAM  USER'S  MANUAL 


A FORTRAN  computer  program  has  been  coded  which  determines  the 
internal  loads,  deflections,  and  relative  stiffness  of  a laminated, 
fixed  ended  beam  with  one  concentrated  load  at  the  center.  The 
program  is  based  on  the  equations  and  procedures  embodied  in  the  pre- 
ceding sections  of  this  report.  A copy  of  the  source  coding  is  con- 
tained in  Appendix  B.  Appendix  C lists  the  intermediate  matrices 
which  were  run  for  Case  A.  The  output  for  four  illustrative  problems 
and  one  test  beam  is  contained  in  Appendix  D.  The  program  was 
written  in  CDC  FORTRAN  Extended  Version  4-0  for  utilization  of  a 
CDC6600  computer. 

Figure  16  shows  a complete  card  deck  ready  to  load  for  a particular 
computer  installation.  Certain  obvious  changes  to  the  information 
shown  on  the  cards,  such  as  time  charge  numbers,  programmer's  name, 
and  program  identification,  would  have  to  be  made  to  suit  the  require- 
ments of  anotner  installation.  A review  of  the  source  coding  shown 
in  Appendix  B may  reveal  other  minor  changes  that  might  be  needed  at 
another  installation. 

The  complete  card  deck  consists  of  essentially  five  categories: 
a job  card,  two  control  cards,  a source  deck,  input  data  cards,  and 
separating  cards,  called  7,  8,  9 Cards  and  6,  7,  8,  9 Cards.  Figure 
17  shows  the  relationship  of  these  components  to  the  complete  deck. 

The  job  card,  control  cards,  and  separating  cards  may  be  obtained 
by  entering  the  appropriate  data  on  the  data  coding  sheet.  Figure  18, 
and  submitting  the  sheet  to  procure  the  key-punched  cards.  The  line 
entries  are  described  as  follows; 

Line  1 - This  line  describes  the  job  card  which  identifies  the 
user  and  lists  other  pertinent  identifying  data. 
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Figure  16.  Complete  card  deck. 


MAIN  PROGRAM 
AND 

SUB-ROUTINES 


Figure  17.  Card  deck  flow  diaqram 


Figure  18.  Sample  coding  sheet  for  control  cards  and  separating  cards 


Line  2 - This  line  refers  to  the  FTN  control  card  and  provides 
access  to  library  subroutines  and  compilation  of  the 
FORTRAN  source  deck. 

Line  3 - This  line  refers  to  an  LGO  control  card,  which  is  a 
load  and  go  card.  One  LGO  card  must  be  inserted  for 
every  set  of  data  cards. 

Line  4 - This  line  refers  to  a card  that  separates  the  source 
deck  from  the  data  cards  and  control  cards.  A minimum 
of  two  cards  is  required  and  an  additional  7,  8,  9 
Card  must  be  added  for  each  additional  set  of  data 
cards. 

Line  5 - This  line  refers  to  the  last  card  in  the  deck,  a 6,  7, 

8,  9 Card. 

The  source  deck  is  a set  of  approximately  1,700  cards  and  may  be 
compiled  from  the  main  program  and  subroutines  listed  in  Appendix  B. 
This  information  may  be  entered  on  data  coding  sheets  and  submitted  to 
key-punch  to  procure  a source  deck. 

The  input  data  card  deck  is  a set  of  six  cards  and  may  be  obtained 
by  entering  the  appropriate  information  on  the  data  coding  sheet. 

Figure  19.  The  data  sheet  may  be  submitted  to  procure  the  key-punched 
cards  shown  in  Figures  20  and  21  for  illustrative  Case  A,  Table  1. 

The  program  statements  are  shown  for  reference.  A set  of  input  data 
cards  must  be  included  for  each  problem  to  be  solved.  The  line  entries 
shown  in  Figure  19  for  Case  A are  described  as  follows: 

Line  1 - Columns  1 through  8 show  the  number  of  spanwise  incre- 
ments to  the  center  line  of  the  beam. 


57 


Figure  19.  Sample  input  data  coding  sheet  (Case  A). 


I 


DATA  CARD  1.  INPUT  NUMBER  OF  SPANWISE  INCREMENTS  TO  CENTER  LINE  OF 
BEAM,  FORCE  (POUNDS),  LENGTH  TO  CENTER  LINE  OF  BEAM 
(INCHES).  WIDTH  (INCHES). 


READ  (5,10)NX,F,L,BB  \ 
10  FORMAT  (I8.3F8.0)  / 


PROGRAM  STATEMENTS  (REFERENCE  ONLY) 


|ll(OtO'0000000'OOOOOOU~'OIOQOBO~OOODODOQOOOOO(IDDOOOOOOODOOOODOOOOI)OOODD000060GOO 

I I I * I I • I I 

inniZiiiiii'Minuimiiii^nniimMiiiimiiMmimimMinniiinuiM 


nnimni 


>21712  inninij 


DATA  CARD  2.  INPUT  THICKNESS  (INCHES)  OF  EACH  LAYER. 


20^F0RMAT°(9F8^0)’'^  " STATEMENTS  (REFERENCE  ONLY) 


0.S5  O.i^ 


0.35  0.  10  0.35  0.  10  0.35 

a : j ■ c 


• 1)0  gg  gggo  ogggooi  oggggog  gg  goog  ggggogg  gg  gggg~ggggogg''tD  ggog'ggggggggggg 

g f I I t t t I I • »f. » 'I  '»•*»!»'  0 ."y  !•  n ft  » P P >'  !•  0 B > •■gii'«g|*»«gi«'«H»«'Uu»io*»ii*g»g 

1 11 11 1 : 11  n t M 11  n n M I n n It  I ^ n n 1 1 n I n n M ; 11 1 1 n n n It  n 1 : 11 1 1 11 1 1 II  t n II  n 


!::nnn'  7 7 2 ]22i  2U 


U 2 2 7 2 ' 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 222  U 


DATA  CARD  3.  INPUT  YOUNG'S  MODULUS  FOR  EACH  LAYER. 


20^F0RMAT°(9F8^0)’'^  " STATEMENTS  (REFERENCE  ONLY) 


/Tooooooo  183.  oi'goooo.  0 133.0340000.0  13s.  Oj4U'iOU.  0 133.  U34000U.  0 

c a a a 3 u a s 


I ■ iiggogoggiQo  g ggggggocg'‘g"ggggggo'gg  B BOCBggB  gg  “g_gBgBgBgg 
; n 1 1 1 1 1 1 1 1 : 1 n 1 n I n n I It  I I n n 1 1 It  1 1 1 1 n c n 1 1 m 1 1 1 1 II 1 1 ;;  II 1 n I n n m n M I n 


1222222222 2221 


>22222222222 22221 


Figure  20.  Sample  data  cards, 
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DATA  CARD  4.  INPUT  SHEAR  MODULUS  FOR  EACH  LAYER. 

READ  (5,20)  (G(J),J  = 1 ,9)  \ PROGRAM  STATEMENTS 
20  FORMAT  (9F8.0)  f (REFERENCE  ONLY) 


1000000. 

n 


73,010^000.0 
P.  c 


73.  0101000.  0 

a C 


73. 01 OluOO.  0 

B B 


73.  P103000.  0 
0 


I" ■"l0»001G0'O  0 G’llllBOBG'I'O O'tOOOOOO'O'J ^D'OCOOIlCfO  G B'lOOGtlOl 

I I I « I i I « I • II  nil  K 'I  ii  i»  'I  .1  j»  »i  n ■> ;«  n m>  I*  N ■ >1  u n »*  n • i • n • «i  «» «i  •'««»•>•  «i  t»  n « im  it  «•'«»•<•  i v u k « • t «i  ii  ••  m n n m ^ n « ^ • 

tniiiniiminriinnnnnin'iiiiMiiiiiinriiiiniiiiuiucmHiMiiiiiii 


( ! ? } 2 2 ? n ! m ! 2 2 2 2 


DATA  CARD  5.  INPUT  DISTANCE  FROM  CENTER  LINE  OF  BEAM  TO  X INCREMENT. 

(NUMBER  OF  VALUES  MUST  EQUAL  NX  ON  CARD  1 AND  LAST  VALUE 
MUST  BE  EQUAL  TO  L ON  CARD  1.)  X(I)  = 0 


READ  (5,20)  (X(J),J  = 2,  NXl))  PROGRAM  STATEMENTS 
20  FORMAT  (9F8.0)  / (REFERENCE  ONLY) 

TTo  TTs  sTo  O O sTs  JTo  TTs  \ 

li  -1  tt  U M j •*  I 


iMio"iooooooflogoo3oooooooooo8~coooDOoooooooo980ooooooaoooDOo~ooooocoooooflosiB 

t I i t I I I I I nriMniiiiiiniiiutiHBii  nn<««iiiian«ti«ti«iiuuHUiiii**4ttii»;uiia«iitn«i'iifiMq«n«i«« 

II  M n M 1 1 1 1 r 1 1 II 1 1 1 ’ 1 1 1 1 1 1 11 1 n I n 1 1 1 1 n n 11 1 1 1 1 1 n 1 1 1 1 1 1 M I n M 1 n M 1 1 1 1 M n 1 


DATA  CARD  6.  (CONTINUATION  OF  CARD  5.) 


5.0 

C 


ll9DOID~fll1l1fiOOOOOOD900000B00800DBIDOB008tlOOflODOtlOQIOOODOOOIMIOOOOOIIDOODOOOOOOOOOO 

I I III  1 I 1 M I 11  1 I I 11  1 I I I 1 I I I I I II  M I I I I I I " I I I I I I I I """""""""""""""  " 

L I J I » I J I J 


! 2 2 2 2 2 2 2 2 2 2 2 1ZJ^ 


~-<m.2  2 2 2 2 2 2 2 2 2 2 2 i 


Figure  21.  Sample  data  cards. 
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Line  1 - Columns  9 through  16  are  the  load  in  pounds. 

Line  1 - Columns  17  through  24  are  the  length  to  the  center  line 
of  the  beam  in  inches. 

Line  1 - Columns  25  through  32  are  the  width  of  the  beam  in 
inches. 

Line  2 - Columns  1 through  8 are  the  thickness  of  the  first  ply 
in  inches  and  the  thickness  of  each  of  the  remaining 
eight  plies  is  shown  in  Columns  9 through  72. 

Line  3 - Columns  1 through  8 are  the  modulus  of  elasticity  (PSI) 
for  Ply  Number  1 and  the  modulfs  of  elasticity  for  Ply 
Numbers  2 through  8 are  shown  in  Columns  9 through  72. 

Line  4 - Columns  1 through  8 are  the  shear  modulus  (PSI)  for 
Ply  Number  1 and  the  shear  modulus  for  Ply  Numbers  2 
through  8 are  shown  in  Columns  9 through  72. 

Lines  5 and  6 - Columns  1 through  8 of  Line  5 are  the  distance 
from  the  center  of  the  beam  t X"  increment  in  inches. 
The  number  of  entries  on  Lines  o and  6 correspond  to 
the  number  of  increments  on  Line  1,  Columns  1 through  8. 
There  are  ten  entries  on  Lines  5 and  6. 

Lines  7 through  12  - These  entries  are  the  appropriate  computer 
symbols  corresponding  to  the  entries  on  Lines  1 through 
6 above  and  are  for  reference  only. 

Data  cards  must  meet  the  following  requirements: 

**  Each  card  field  is  eight  (8)  card  columns  wide. 
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® There  cannot  be  any  zeros  or  blank  fields  imbedded 
in  the  data  (i.e.,  If  a layer  is  to  be  dropped  from 
the  problem,  its  thickness  and  other  parameters 
must  be  input  to  the  problem  as  relatively  small 
numbers  to  cause  its  effective  elimination.)* 

" Numbers  must  be  right  justified  in  their  fields. 

® A trailing  decimal  point  may  be  dropped. 
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SECTION  IV 
CONCLUSIONS 


The  analysis  and  computer  program  presented  in  this  report  were 
developed  to  provide  a tool  for  evaluating  laminated  combinations  of 
transparent  materials  during  initial  design  of  windshields  and  windows. 
Historically,  formulas  and  theories  from  engineering  handbooks  have 
been  used  for  this  purpose,  but  these  methods  are  highly  approximate. 

The  laminated  beam  analyzed  in  the  present  approach  is  considered 
to  represent  a strip  cut  from  a typical  aircraft  windshield  transparency. 
The  primary  characteristic  of  such  a laminated  beam  is  the  relative  ease 
of  structural  plies  to  slide  past  each  other  because  of  the  softness  of 
interlayer  materials.  Equilibrium  and  compatibility  equations  are 
written  for  the  beam,  based  upon  assumptions  that  preserve  the  important 
features  of  laminated  beam  behavior.  The  resulting  set  of  differential 
equations  is  solved  exactly. 

The  computer  code  presented  can  expedite  the  application  of  the 
theory  to  practical  problems.  Utilizing  this  program,  the  deflections 
for  a series  of  typical  laminated  beams  were  examined,  and  were  compared 
to  the  results  of  a proven  finite  element  computer  program.  The  two 
sets  of  computed  values  agreed  within  one  percent.  Data  obtained  <^rom 
another  beam,  tested  as  noted  in  Reference  1,  were  compared  to  calcu- 
lated results  provided  by  the  present  method.  The  outer  ply  of  the 
test  beam  was  "free  floating"  at  the  ends,  rather  than  fixed,  as 
assumed  in  the  analysis.  To  simulate  this  condition,  a negligible  value 
of  Young's  modulus  was  assigned  to  the  outer  ply.  The  calculated  deflec- 
tion was  3.5  percent  greater  than  the  measured  value.  This  agreement  is 
considered  reasonably  close,  in  view  of  the  approximate  method  used  to 
account  for  the  "free  floating"  outer  ply. 

The  applications  of  the  computer  program  described  in  Section  II 
under  the  heading  "Illustrative  Examples"  were  run  on  the  CDC6600 
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computer,  employing  21,000  decimal  words  of  core.  Each  case  required 
0.879  seconds  CPU  time  and  7.588  seconds  I/O  time.  The  I/O  time 
Includes  compilation  of  the  entire  program. 
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SECTION  V 
RECOMMENDATIONS 

The  new  computer  program  is  recommended  as  a tool  for  use  in  the 
early  stages  of  aircraft  windshield  design  as  a means  of  screening  lami- 
nate configurations.  As  an  example,  a strip  of  a certain  width  can  be 
cut  out  of  a transparency  being  studied.  This  strip  can  be  flattened 
to  form  a laminated  beam  having  pinned  or  fixed  ends  as  appropriate. 

The  width  of  the  strip  can  be  selected  so  that  the  resulting  beam,  sup- 
ported only  at  the  ends,  roughly  represents  the  transparency  supported 
on  all  edges.  A load  can  be  applied  at  the  center  to  represent  a bird 
impact.  The  number  of  plies,  structural  ply  and  interlayer  materials, 
and  ply  thicknesses  can  be  varied,  so  that  establisheu  stress  and/or 
displacement  constraints  are  met.  The  appropriate  constraints  are 
that  ultimate  stresses  and  maximum  allowable  displacements  are  not 
exceeded.  Weights  per  square  foot  of  transparency  surface  can  then  be 
calculated  and  compared.  The  laminates  showing  the  most  promise  on 
this  basis  can  be  selected  for  further  study. 

Another  use  for  the  present  code  is  to  provide  an  effective  bending 
stiffness  for  a laminate.  This  effective  stiffness  can  be  selected  as 
the  El  value  of  a monolithic  beam  necessary  to  give  the  same  maximum 
deflection  as  a beam  composed  of  the  given  laminate  under  the  same  load. 
The  computer  program  outputs  such  a value.  This  value  can  then  be 
employed  as  a means  of  developing  a finite  element  model  of  the  laminated 
transparency  and  supporting  structure,  in  which  the  laminate  is  replaced 
by  a single  layer,  as  a means  of  reducing  engineering,  computing,  and 
elapsed  time  to  perform  either  a static  or  a dynamic  analysis.  Since 
effective  stiffness  depends  upon  beam  length,  the  user  should  estimate 
the  half  length  of  the  wave  produced  in  the  transparency  when  impact 
occurs.  Such  an  estimate  can  be  made  on  the  basis  of  past  experience, 
including  previous  studies  employing  the  Bird  Impact  Math  Model 
(Reference  6)  or  test  results.  This  half  wave  length  can  then  be  taken 
as  the  length  of  the  simply  supported  beam  for  which  an  effective  El 
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value  can  be  calculated.  An  approach  of  this  kind  is  probably  most  use- 
ful for  calculating  transparency  displacements  and  loads  on  edge  attach- 
ments and  supporting  structure. 

The  feasibility  of  extending  the  method  to  cover  other  loadings  and 
boundary  conditions  should  be  studied.  For  example,  an  option  applicable 
to  an  end  condition  in  which  the  layers  are  prevented  from  sliding  rela- 
tive to  each  other,  as  by  a bolt,  although  the  end  tangent  can  rotate, 
might  be  useful.  A beam  uniformly  loaded  over  part  of  its  length  is  a 
candidate  case.  The  feasibility  of  applying  the  approach  to  dynamic 
problems  also  should  be  studied,  and  the  possibility  of  eliminating  the 
assumption  that  transverse  strains  are  negligible  should  be  considered. 
Eliminating  this  assumption  would  mean  that  the  method  would  provide 
reasonably  reliable  estimates  of  transverse  stresses,  and  perhaps  a 
means  of  forecasting  delaraination.  Any  of  these  improvements  would 
increase  the  realism  of  the  analytical  results.  Consequently,  the 
usefulness  of  the  approach  as  a design  tool  would  be  enhanced. 

If  any  of  these  improvements  were  accomplished,  they  could  be 
verified  by  comparison  with  results  obtained  by  finite  element  analysis. 
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APPENDIX  A 


PROPERTIES  OF  THE  EIGENVALUES 

and  eigenvectors 


The  matric  differential  equation 
in  the  following  form  by  rearranging 


(Equation  42)  can  be  written 
the  rows  of  Y : 
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The  corresponding  characteristic  equation  is 


where 
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(A. 2) 


(A. 3) 

(A. 4) 
(A. 5) 
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Changing  the  signs  of  x and  6|^  leaves  Equations  A. 5 and  A. 6 
unchanged. 

vector  of  Equation  A. 4,  then  -x 
eigenvalue  and  an  eigenvector. 


Therefore  if  x and 
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j are  an  eigenvalue  and  an  eigen- 
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are  also  an 
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SOURCE  CODING  (MAIN  PROGRAM  AND  SUB-ROUTINES) 
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IS  DESTROYED.  15490 
IF  MORE  THAN  30  ITERATIONS  ARE  REQUIRED  TO  DETERMINE  15500 
AN  EIGENVALUE,  THIS  SUBROUTINE  TERMINATES  WITH  IND  15510 
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THE  EIGENVECTORS  ARE  DETERMINED  BY  BACK-SUBSTITUTION  AND  15800 
THEN  TRANSFORMED  TO  THE  EIGENVECTORS  OF  THE  ORIGINAL  15810 
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K = 1 16210 
IF  (EIGVAL(2,K).NE.0.0)  GO  TO  6 16220 
DO  5 J=1  .M  16230 


K = I 16530 

CMOD  = TEMP  16540 

13  CONTINUE  16550 

DENOM  = CMPLX(EIGVEC(1 ,K,J),EIGVEC(2,K,J))  16560 
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